Prac5_London_airbnb_analysis

# loading packages/libraries
library(sf)

## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

library(here)

## here() starts at D:/CASA/GIS/gis_code/Prac5

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(tmap)

## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')

library(tmaptools)
library(grid)
library(leafpop)
library(leaflet)

1. Spatial Join

# reading in the London Boroughs data
Londonborough <- st_read(here("prac5_data", 
                              "statistical-gis-boundaries-london",
                              "ESRI",
                              "London_Borough_Excluding_MHW.shp")) %>%
  select(c(-SUB_2009,-SUB_2006)) %>%
  st_transform(., 27700) # transforming the CRS of this data set

## Reading layer `London_Borough_Excluding_MHW' from data source 
##   `D:\CASA\GIS\gis_code\Prac5\prac5_data\statistical-gis-boundaries-london\ESRI\London_Borough_Excluding_MHW.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 33 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
## Projected CRS: OSGB36 / British National Grid

# reading the OSM data
OSM <- st_read(here("prac5_data", 
                    "greater-london-latest-free.shp", 
                    "gis_osm_pois_a_free_1.shp")) %>%
    st_transform(., 27700) %>%
  #select hotels only
  dplyr::filter(fclass == 'hotel')

## Reading layer `gis_osm_pois_a_free_1' from data source 
##   `D:\CASA\GIS\gis_code\Prac5\prac5_data\greater-london-latest-free.shp\gis_osm_pois_a_free_1.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 48312 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -0.5108706 ymin: 51.28109 xmax: 0.322123 ymax: 51.6926
## Geodetic CRS:  WGS 84

# a  spatial join example for these two datasets
join_example <-  st_join(Londonborough, OSM)

head(join_example)

## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 516362.6 ymin: 159907.4 xmax: 522654 ymax: 172367
## Projected CRS: OSGB36 / British National Grid
##                     NAME  GSS_CODE HECTARES NONLD_AREA ONS_INNER    osm_id code
## 1   Kingston upon Thames E09000021 3726.117          0         F  28055159 2401
## 1.1 Kingston upon Thames E09000021 3726.117          0         F 161631526 2401
## 1.2 Kingston upon Thames E09000021 3726.117          0         F 421483573 2401
## 1.3 Kingston upon Thames E09000021 3726.117          0         F 421506999 2401
## 1.4 Kingston upon Thames E09000021 3726.117          0         F 502772206 2401
## 1.5 Kingston upon Thames E09000021 3726.117          0         F 893668667 2401
##     fclass                                      name
## 1    hotel                               Premier Inn
## 1.1  hotel                  Chessington Safari Hotel
## 1.2  hotel                   Premier Inn Chessington
## 1.3  hotel                  Chessington Azteca Hotel
## 1.4  hotel                              Warren House
## 1.5  hotel DoubleTree by Hilton Kingston Upon Thames
##                           geometry
## 1   MULTIPOLYGON (((516401.6 16...
## 1.1 MULTIPOLYGON (((516401.6 16...
## 1.2 MULTIPOLYGON (((516401.6 16...
## 1.3 MULTIPOLYGON (((516401.6 16...
## 1.4 MULTIPOLYGON (((516401.6 16...
## 1.5 MULTIPOLYGON (((516401.6 16...

2. Static Map

# won't run because disabled by eval=FALSE, and because I have done this above only here for notes purposes as we shift into Static Map
#Load all our data
library(sf)
library(tmap)
library(tmaptools)
library(tidyverse)
library(here)

# read in all the spatial data and 
# reproject it 

OSM <- st_read(here::here("prac5_data",
                          "greater-london-latest-free.shp", 
                          "gis_osm_pois_a_free_1.shp")) %>%
  st_transform(., 27700) %>%
  #select hotels only
  filter(fclass == 'hotel')

##London Borough data is already in 277000
Londonborough <- st_read(here::here("Prac1_data",
                                    "statistical-gis-boundaries-london", 
                                    "ESRI", 
                                    "London_Borough_Excluding_MHW.shp"))%>%
  st_transform(., 27700)

# loading more data
Worldcities <- st_read(here::here("prac5_data", 
                                  "World_Cities", 
                                  "World_Cities.shp")) %>%
  st_transform(., 27700)

## Reading layer `World_Cities' from data source 
##   `D:\CASA\GIS\gis_code\Prac5\prac5_data\World_Cities\World_Cities.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2540 features and 14 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -19609100 ymin: -7321602 xmax: 19950890 ymax: 14476660
## Projected CRS: WGS 84 / Pseudo-Mercator

UK_outline <- st_read(here::here("prac5_data", 
                                 "gadm41_GBR_shp", 
                                 "gadm41_GBR_0.shp")) %>%
  st_transform(., 27700)

## Reading layer `gadm41_GBR_0' from data source 
##   `D:\CASA\GIS\gis_code\Prac5\prac5_data\gadm41_GBR_shp\gadm41_GBR_0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -8.649996 ymin: 49.86531 xmax: 1.764393 ymax: 60.84548
## Geodetic CRS:  WGS 84

Example code to unzip a .gz file

# read in the .csv
# and make it into spatial data

Airbnb <- read_csv(here("prac5_data", 
                        "listings.csv")) %>%
  # longitude is considered x value here, latitude is y
  st_as_sf(., 
           coords = c("longitude", "latitude"), 
           crs = 4326) %>%
    st_transform(., 27700)%>%
    #select entire places that are available all year
    filter(room_type == 'Entire home/apt' & availability_365 =='365')

## Rows: 96182 Columns: 18
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr   (4): name, host_name, neighbourhood, room_type
## dbl  (11): id, host_id, latitude, longitude, price, minimum_nights, number_o...
## lgl   (2): neighbourhood_group, license
## date  (1): last_review
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Lets make a function to join our data. It will aggregate counts of the number of times a borough appears

# make a function for the join
# functions are covered in practical 7
# but see if you can work out what is going on
# hint all you have to do is replace data1 and data2
# with the data you want to use

Joinfun <- function(data1, data2){

output<- data1%>%
  st_join(data2,.) %>% # notice the second data set is on the left
  add_count(GSS_CODE, name="hotels_in_borough") # function creating counts of GSS_CODE and storing in column "hotels_in_borough" 

  return(output)
}

Use the new function to join our data sets

# use the function for hotels
Hotels <- Joinfun(OSM, Londonborough)

# then for airbnb
# this is incorrect - change to airbnb2 to look at result
Airbnb <- Joinfun(Airbnb, Londonborough)

Worldcities2 <- Worldcities %>%
  filter(CNTRY_NAME=='United Kingdom'&
           Worldcities$CITY_NAME=='Birmingham'|
           Worldcities$CITY_NAME=='London'|
           Worldcities$CITY_NAME=='Edinburgh')

newbb <- c(xmin=-296000, ymin=5408, xmax=655696, ymax=1000000) # setting a bounding box
  
UK_outlinecrop <- UK_outline$geometry %>%
  st_crop(., newbb) # doing a spatial crop using the new bounding box

Hotels <- Hotels %>%
  #at the moment each hotel is a row for the borough
  #we just want one row that has number of airbnbs
  group_by(., GSS_CODE, NAME) %>%
  summarise(`Accomodation count` = unique(hotels_in_borough))

## `summarise()` has grouped output by 'GSS_CODE'. You can override using the
## `.groups` argument.

Airbnb <- Airbnb %>%
  group_by(., GSS_CODE, NAME)%>%
  summarise(`Accomodation count` = unique(hotels_in_borough))

## `summarise()` has grouped output by 'GSS_CODE'. You can override using the
## `.groups` argument.

Hotels %>%
  filter(NAME=="Sutton")

## Simple feature collection with 1 feature and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 522223.3 ymin: 159658.1 xmax: 531245.8 ymax: 167622.5
## Projected CRS: OSGB36 / British National Grid
## # A tibble: 1 × 4
## # Groups:   GSS_CODE [1]
##   GSS_CODE  NAME   `Accomodation count`                                 geometry
## * <chr>     <chr>                 <int>                       <MULTIPOLYGON [m]>
## 1 E09000029 Sutton                    1 (((528552.3 159658.1, 528399.7 159928.8…

Making the map

tmap_mode("plot")

## tmap mode set to plotting

# set the breaks for our mapped data
breaks = c(0, 5, 12, 26, 57, 286) 

# plot each map
tm1 <- tm_shape(Hotels) + 
  tm_polygons("Accomodation count", 
              breaks=breaks,
              palette="PuBu")+
  tm_legend(show=FALSE)+
  tm_layout(frame=FALSE)+
  tm_credits("(a)", position=c(0,0.85), size=1.5)

tm2 <- tm_shape(Airbnb) + 
  tm_polygons("Accomodation count",
              breaks=breaks, 
              palette="PuBu") + 
  tm_legend(show=FALSE)+
  tm_layout(frame=FALSE)+
  tm_credits("(b)", position=c(0,0.85), size=1.5)

tm3 <- tm_shape(UK_outlinecrop)+ 
  tm_polygons(col="darkslategray1")+
  tm_layout(frame=FALSE)+
  tm_shape(Worldcities2) +
  tm_symbols(col = "red", scale = .5)+
  tm_text("CITY_NAME", xmod=-1, ymod=-0.5)

legend <- tm_shape(Hotels) +
    tm_polygons("Accomodation count",
                breaks=breaks,
                palette="PuBu") +
    tm_scale_bar(position=c(0.2,0.04), text.size=0.6)+
    tm_compass(north=0, position=c(0.65,0.6))+
    tm_layout(legend.only = TRUE, legend.position=c(0.2,0.25),asp=0.1)+
    tm_credits("(c) OpenStreetMap contrbutors and Air b n b", position=c(0.0,0.0))
  
t=tmap_arrange(tm1, tm2, tm3, legend, ncol=2)

t

## Warning: Values have found that are higher than the highest break
## Warning: Values have found that are higher than the highest break

Arranging maps using the grid package

#library(grid)
# erases the current device or moves to a new page 
# probably not needed but makes sure you are plotting on a new page.
grid.newpage()

pushViewport(viewport(layout=grid.layout(2,2)))
print(tm1, vp=viewport(layout.pos.col=1, layout.pos.row=1, height=5))

## Warning: Values have found that are higher than the highest break

print(tm2, vp=viewport(layout.pos.col=2, layout.pos.row=1, height=5))
print(tm3, vp=viewport(layout.pos.col=1, layout.pos.row=2, height=5))
print(legend, vp=viewport(layout.pos.col=2, layout.pos.row=2, height=5))

## Warning: Values have found that are higher than the highest break

Inset Map

# calculating the bounding box for the London Airbnb dataset
Londonbb <- st_bbox(Airbnb,
                    crs = st_crs(Airbnb))%>%
  #we need this to convert it into a class of sf
  # otherwise if our bb won't have a class it will just be x and y coordinates for the box
  # this makes it into a polygon
  st_as_sfc()


main <- tm_shape(Airbnb, bbbox = Londonbb) + 
  tm_polygons("Accomodation count",
              breaks=breaks, 
              palette="PuBu")+
  tm_scale_bar(position = c("left", "bottom"), text.size = .75)+
  tm_layout(legend.position = c("right","top"), 
            legend.text.size=.75, 
            legend.title.size = 1.1,
            frame=FALSE)+
  tm_credits("(c) OpenStreetMap contrbutors and Air b n b", position=c(0.0,0.0))+
  #tm_text(text = "NAME", size = .5, along.lines =T, remove.overlap=T,  auto.placement=F)+
  tm_compass(type = "8star", position = c(0.06, 0.1)) +

  #bottom left top right <- all round margin
  tm_layout(inner.margin=c(0.02,0.02,0.02,0.2))


inset = tm_shape(UK_outlinecrop) + 
  tm_polygons() +
  tm_shape(Londonbb)+ 
  tm_borders(col = "grey40", lwd = 3)+ ##lwd==line width
    tm_layout(frame=FALSE,
            bg.color = "transparent")+
  tm_shape(Worldcities2) +
  tm_symbols(col = "red", scale = .5)+
  tm_text("CITY_NAME", xmod=-1.5, ymod=-0.5) # positioning of the text with xmod and ymod, -ve means left and down respectively


#library(grid)
main
print(inset, vp = viewport(0.86, 0.29, width = 0.5, height = 0.55))

Exporting

tmap_save(t, 'hotelsandairbnbR.png')

## Warning: Values have found that are higher than the highest break
## Warning: Values have found that are higher than the highest break

## Map saved to D:\CASA\GIS\gis_code\Prac5\hotelsandairbnbR.png

## Resolution: 2100 by 2100 pixels

## Size: 7 by 7 inches (300 dpi)

#library(grid)
tmap_save(tm = main,insets_tm = inset, insets_vp = viewport(0.86, 0.29, width=.5, height=.55), filename="test.pdf", dpi=600)

## Map saved to D:\CASA\GIS\gis_code\Prac5\test.pdf

## Size: 8.833333 by 5.541667 inches

Interactive maps

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(Airbnb) + 
  tm_polygons("Accomodation count", breaks=breaks)

Advanced Interactive map

# library for pop up boxes
#library(leafpop)
#library(leaflet)

#join data
Joined <- Airbnb%>%
  st_join(., Hotels, join = st_equals)%>%
  dplyr::select(GSS_CODE.x, NAME.x, `Accomodation count.x`, `Accomodation count.y`)%>%
  dplyr::rename(`GSS code` =`GSS_CODE.x`,
                `Borough` = `NAME.x`,
                `Airbnb count` = `Accomodation count.x`,
                `Hotel count`= `Accomodation count.y`)%>%
  st_transform(., 4326)
  
  
#remove the geometry for our pop up boxes to avoid
popupairbnb <-Joined %>%
  st_drop_geometry()%>%
  dplyr::select(`Airbnb count`, Borough)%>%
  popupTable()

## Adding missing grouping variables: `GSS code`

popuphotel <-Joined %>%
  st_drop_geometry()%>%
  dplyr::select(`Hotel count`, Borough)%>%
  popupTable()

## Adding missing grouping variables: `GSS code`

tmap_mode("view")

## tmap mode set to interactive viewing

# set the colour palettes using our previously defined breaks


pal1 <- Joined %>%
  colorBin(palette = "YlOrRd", domain=.$`Airbnb count`, bins=breaks)

pal1 <-colorBin(palette = "YlOrRd", domain=Joined$`Airbnb count`, bins=breaks)

pal2 <- Joined %>%
  colorBin(palette = "YlOrRd", domain=.$`Hotel count`, bins=breaks)


map<- leaflet(Joined) %>%

  #add our polygons, linking to the tables we just made
  addPolygons(color="white", 
              weight = 2,
              opacity = 1,
              dashArray = "3",
              popup = popupairbnb,
              fillOpacity = 0.7,
              fillColor = ~pal2(`Airbnb count`),
              group = "Airbnb")%>%
  
  addPolygons(fillColor = ~pal2(`Hotel count`), 
              weight = 2,
              opacity = 1,
              color = "white",
              dashArray = "3",
              popup = popupairbnb,
              fillOpacity = 0.7,group = "Hotels")%>%
  
  #add basemaps
  addTiles(group = "OSM (default)") %>%
  addProviderTiles(providers$Stadia.StamenToner, group = "Toner") %>%
  addProviderTiles(providers$Stadia.StamenTonerLite, group = "Toner Lite") %>%
  addProviderTiles(providers$CartoDB.Positron, group = "CartoDB")%>%
  
  # add a legend
  addLegend(pal = pal2, values = ~`Hotel count`, group = c("Airbnb","Hotel"), 
            position ="bottomleft", title = "Accomodation count") %>%
  # specify layers control
  addLayersControl(
    baseGroups = c("OSM (default)", "Toner", "Toner Lite", "CartoDB"),
    overlayGroups = c("Airbnb", "Hotels"),
    options = layersControlOptions(collapsed = FALSE)
  )

## Warning in pal2(`Hotel count`): Some values were outside the color scale and
## will be treated as NA

# plot the map
map

# this is what was explained in Prac6 section: 6.5.4
# attempting a spatial join, by default, it uses st_intersects which is problematic. See below
all_accomodation <- st_join(Hotels,Airbnb)

head(all_accomodation)

## Simple feature collection with 6 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 530967.7 ymin: 180510.3 xmax: 533839.6 ymax: 182206.1
## Projected CRS: OSGB36 / British National Grid
## # A tibble: 6 × 7
## # Groups:   GSS_CODE.x [1]
##   GSS_CODE.x NAME.x  `Accomodation count.x`                  geometry GSS_CODE.y
##   <chr>      <chr>                    <int>        <MULTIPOLYGON [m]> <chr>     
## 1 E09000001  City o…                     23 (((531145.1 180782.1, 53… E09000001 
## 2 E09000001  City o…                     23 (((531145.1 180782.1, 53… E09000007 
## 3 E09000001  City o…                     23 (((531145.1 180782.1, 53… E09000012 
## 4 E09000001  City o…                     23 (((531145.1 180782.1, 53… E09000019 
## 5 E09000001  City o…                     23 (((531145.1 180782.1, 53… E09000030 
## 6 E09000001  City o…                     23 (((531145.1 180782.1, 53… E09000033 
## # ℹ 2 more variables: NAME.y <chr>, `Accomodation count.y` <int>

# lets correct this by changing the argument to st_equals
all_accomodation2 <- st_join(Hotels, Airbnb, join = st_equals)

head(all_accomodation2)

## Simple feature collection with 6 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 515484.9 ymin: 156480.9 xmax: 554089.2 ymax: 198355.2
## Projected CRS: OSGB36 / British National Grid
## # A tibble: 6 × 7
## # Groups:   GSS_CODE.x [6]
##   GSS_CODE.x NAME.x  `Accomodation count.x`                  geometry GSS_CODE.y
##   <chr>      <chr>                    <int>        <MULTIPOLYGON [m]> <chr>     
## 1 E09000001  City o…                     23 (((531145.1 180782.1, 53… E09000001 
## 2 E09000002  Barkin…                      9 (((543905.4 183199.1, 54… E09000002 
## 3 E09000003  Barnet                      14 (((524579.9 198355.2, 52… E09000003 
## 4 E09000004  Bexley                       6 (((547226.2 181299.3, 54… E09000004 
## 5 E09000005  Brent                       13 (((525201 182512.6, 5251… E09000005 
## 6 E09000006  Bromley                      4 (((540373.6 157530.4, 54… E09000006 
## # ℹ 2 more variables: NAME.y <chr>, `Accomodation count.y` <int>